function ic_fsw_script(n, L, ampl_z, ampl_psi, k0, kc, tmax, cfl, h, wskip, apow, av_coef, grav, fname)
%function ic_fsw_script(n, L, ampl, k0, kc, tmax, cfl, h, wskip, apow, av_coef, grav, fname) generates IC in the form of gaussian 
%on n gridpoints on the interval [-L/2,L/2] with amplitude ampl, centered at k0 and width kc
%and writes it to file fname. Will evolve until tmax with dt = cfl*L/n at water depth h.
%Profiles written every wskip seconds. Artificial viscosity of power apow with coefficient av_coef
%Gravity coeff set to grav.
%
%
%
%
  format long e
  IC = ampl_psi; IC2 = ampl_z;
  dx = L/n; x = (0.5:n-0.5)*dx - L/2; 
  c = 2*pi/L; k = c*(-n/2:n/2-1); 


  z = zeros(1,n);
  A = IC.*exp(-((k-k0)/kc).^2);
  A2 = IC2.*exp(-(k/kc).^2);

  z = A2.*exp(2i*pi*rand(1,n));
  
  psi = A.*exp(2i*pi*rand(1,n));


  zr2 = ifft(ifftshift(z))*sqrt(2*c)*n/L;
  psi2 = real(ifft(ifftshift(psi))*sqrt(2*c)*n/L);


  %psi2 = 0.4*cos(c*x);
  %zr2 = 0.*exp(i*c*x);

  max(abs(psi2))
  fps = fftshift(fft(psi2));
  k_energy = sum(0.5*(k.*tanh(k*h)).*(conj(fps).*fps))/n^2


  dz = 1+ifft(ifftshift(1i*k.*fftshift(fft(zr2))));
  max(imag(zr2))
  p_energy = sum(0.5*grav*(imag(zr2).^2 ).*(real(dz)))*dx

  figure(2)
  semilogy(k,[A; A2],0.66*max(k),1e-4,'+');
  axis([min(k) max(k) 1e-15 IC])
  figure(2)
  subplot(2,1,1)
  plot(x,[real(psi2); imag(psi2)]);
  subplot(2,1,2)
  plot(x,[real(zr2); imag(zr2)]);

  fh = fopen(fname,'w');
  fprintf(fh, '## Input File for free surface wave run\n#\n#gridpoints\t%d\n', n);
  fprintf(fh, '#max_time\t%f\n#cfl_num\t%f\n#threads\t%d\n#\n', tmax, cfl, 2);
  fprintf(fh, '##Info to generate IC\n#\n');
  fprintf(fh, '#box_length\t%f\n#water_depth\t%f\n', L, h);
  fprintf(fh, '#ic_type\t%d\n#height_psi\t%f\n#height_z\t%f\n', 1, 1, 1);
  fprintf(fh, '#wavenum_psi\t%e\n#wavenum_z\t%e\n', 3, 5 );
  fprintf(fh, '#height_z2\t%e\n#wavenum_z2\t%e\n', 5e-1, 2e-1);
  fprintf(fh, '#write_skip\t%e\n#cutoff_knum\t%d\n', wskip, 128);
  fprintf(fh, '#art_viscos\t%e\n#grav_accel\t%e\n', av_coef, grav);
  fprintf(fh, '#k_power\t%d\n#epsilon\t%e\n\n', apow, 0);
  fprintf(fh, '');
  fprintf(fh, '');
  fprintf(fh, '## IC generated by ./ic_fsw_script.m\n');
  fprintf(fh, '#1. x\t\t2. Psi_re\t\t3.Psi_im\t\t4. z_re\t\t5. z_im\n\n');
  


  for j = 1:n
    fprintf(fh, '%e\t%.15e\t%.15e\t%.15e\t%.15e\n',x(j), real(psi2(j)),imag(psi2(j)), real(zr2(j)),imag(zr2(j)));
  end 
  fclose(fh);
end
